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Abstract 

a-stable distributions are utilised as models for heavy-tailed noise in many 
areas of statistics, finance and signal processing engineering. However, in gen- 
eral, neither univariate nor multivariate a-stable models admit closed form 
densities which can be evaluated pointwise. This complicates the inferential 
procedure. As a result, a-stable models are practically limited to the univari- 
ate setting under the Bayesian paradigm, and to bivariate models under the 
classical framework. In this article we develop a novel Bayesian approach to 
modelling univariate and multivariate a-stable distributions based on recent 
advances in "likelihood-free" inference. We present an evaluation of the per- 
formance of this procedure in 1, 2 and 3 dimensions, and provide an analysis 
of real daily currency exchange rate data. The proposed approach provides a 
feasible inferential methodology at a moderate computational cost. 

Keywords: a-stable distributions; Approximate Bayesian computation; Likelihood- 
free inference; Sequential Monte Carlo samplers. 
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1 Introduction 



Models constructed with a-stable distributions possess several useful properties, in- 
cluding infinite variance, skewness and heavy tails (IZolotarev 1986t lAlder et al. 1998} 
Samorodnitsky and Taqqu 1994 INolan 20071) . a-stable distributions provide no gen- 



eral analytic expressions for the density, median, mode or entropy, but are uniquely 
specified by their characteristic function, which has several parameterizations. Con- 
sidered as generalizations of the Gaussian distribution, they are defined as the class 
of location-scale distributions which are closed under convolutions, a-stable distribu- 
tions have found application in many areas of statistics, finance and signal processing 
engineering as models for impulsive, heavy tailed noise processes ( IMandelbrot 1960 



IFama 19651 [Fama and Roll 19681 iNiEas and Shao 19951 IGodsill"2000llMelchiori 2006|1 



The univariate a-stable distribution is typically specified by four parameters: a G 
(0, 2] determining the rate of tail decay; /3 G [—1,1] determining the degree and sign of 
asymmetry (skewness); 7 > the scale (under some parameterizations); and (5 G M the 



location ( Levy 1924 ). The parameter a. is termed the characteristic exponent, with 
small and large a implying heavy and light tails respectively. Gaussian (a = 2, /3 = 0) 
and Cauchy (a = l,/3 = 0) distributions provide the only analytically tractable sub- 
members of this family. In general, as a-stable models admit no closed form expression 
for the density which can be evaluated pointwise (excepting Gaussian and Cauchy 
members), inference typically proceeds via the characteristic function. 

This paper is concerned with constructing both univariate and multivariate Bayesian 
models in which the likelihood model is from the class of a-stable distributions. This 
is known to be a difficult problem. Existing methods for Bayesian a-stable mod- 
els are limited to the univariate setting (IBuckle 1995t IGodsill 1999t IGodsill 2000t 
ILombardi 20071 ICasarin^OOil ISalas-Gonzalez et al. 20061) . 
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Inferential procedures for a-stable models may be classified as auxiliary variable 
methods, inversion plus series expansion approaches and density estimation methods. 
The auxiliary variable Gibbs sampler fIBuckle 1995^ increases the dimension of the 
parameter space from 4 (a, /3, 7 and 6) to n+4, where n is the number of observations. 
As strong correlations between parameters and large sample sizes are common in the 
a-stable setting, this results in a slowly mixing Markov chain since Gibbs moves 
are limited to moving parallel to the axes (e.g. INeal 20031) . Other Markov chain 
Monte Carlo (MCMC) samplers (IDuMouchel 19751 ILombardi 20071) adopt inversion 
techniques for numerical integration of the characteristic function, employing inverse 



Fourier transforms combined with a series expansion (Bergstrom 1953) to accurately 
estimate distributional tails. This is performed at each iteration of the Markov chain 
to evaluate the likelihood, and is accordingly computationally intensive. In addition 
the quality of the resulting approximation is sensitive to the spacing of the fast Fourier 
transform grid and the point at which the series expansion begins (ILombardi 2007^ . 

Univariate density estimation methods include integral representations (IZolotarev 19861) . 
parametric mixtures (INolan 1997} IMcCuUoch 1998|) and numerical estimation through 



splines and series expansions ( INolan et al. 200"T| INolan 20081) . McCuUoch (1998) ap- 
proximates symmetric stable distributions using a mixture of Gaussian and Cauchy 



densities. Doganoglu and Mittnik (1998) and Mittnik and Rachev (1991) approxi- 



mate the a-stable density through spline polynomials, and Kuruoglu et al. (1997) 



via a mixture of Gaussian distributions. Parameter estimation has been performed 
by an expectation-maximization (EM) algorithm (ILombardi and Godsill 2006P and 
by method of (partial) moments (IPress 1972t IWeron 2006p . Implemented within an 
MCMC sampler, such density estimation methods would be highly computational. 

None of the above methods easily generalize to the multivariate setting. It is 
currently only practical to numerically evaluate two dimensional a-stable densities via 
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inversion of the characteristic function. Here the required computation is a function 
of a and the number and spread of masses in the discrete spectral representation 
( INolan et aL 200"T|: INolan 1997p . Beyond two dimensions this procedure becomes 
untenably slow with limited accuracy. 

In this article we develop practical Bayesian inferential methods to fit univariate 
and multivariate a-stable models. To the best of our knowledge, no practical Bayesian 
methods have been developed for the multivariate model as the required computa- 
tional complexity increases dramatically with model dimension. The same is true of 
classical methods beyond two dimensions. Our approach is based on recent develop- 
ments in "likelihood-free" inference, which permits approximate posterior simulation 
for Bayesian models without the need to explicitly evaluate the likelihood. 

In Section [2] we briefly introduce likelihood- free inference and the sampling frame- 
work used in this article. Section [3] presents the Bayesian a-stable model, with a par- 
ticular focus on summary statistic specification, a critical component of likelihood-free 
inference. We provide an evaluation of the performance of the proposed methodology 
in Section m based on controlled simulation studies in 1, 2 and 3 dimensions. Fi- 
nally, in Section [5] we demonstrate an analysis of real daily currency data under both 
univariate and multivariate settings, and provide comparisons with existing methods. 
We conclude with a discussion. 

2 Likelihood-free models 

Computational procedures to simulate from posterior distributions, T^{Q\y) oc 7r(y|6')7r(6'), 
of parameters ^ G 6 given observed data y G A", are well established (e.g. IBrooks et al. 20T0l) . 
However when pointwise evaluation of the likelihood function 7r(y|6') is computa- 
tionally prohibitive or intractable, alternative procedures are required. Likelihood- 
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free methods (also known as approximate Bayesian computation) permit simulation 
from an approximate posterior model while circumventing explicit evaluation of the 
likelihood function fITavare et al. 1997} IBeaumont et al. 2002t Marjoram et al. 2003 
ISisson et al. 20071 [Ratmann et al. 2009|) . 



Assuming data simulation x ~ 7r(x|6') under the model given 6 is easily obtainable, 
likelihood-free methods embed the posterior 7r{y\6) within an augmented model 

TTLFiO, x\y) cc 7T,iy\x, e)7r{x\9)n{e), (2.1) 

where x ~ 7r(a;|6'), x G X, is an auxiliary parameter on the same space as the 
observed data y. The function TT^{y\x,6) is typically a standard smoothing kernel 
(e.g. IBlum 20091) with scale parameter e, which weights the intractable posterior 
with high values in regions when the observed data y and auxiliary data x are simi- 
lar. For example, uniform kernels are commonplace in likelihood-free models (e.g. 
Marjoram et al. 2003 ISisson et al. 20071) . although alternatives such as Epanech- 



nikov (IBeaumont et al. 2002^ and Gaussian kernels (IPeters et al. 20081) provide im- 
proved efficiency. The resulting approximation to the true posterior target distribu- 
tion 

7rLF{0\y)cx [ 7rM^,e)n{x\e)n{e)dx = 7T{e)E^^^ie)[rre{y\x,e)] (2.2) 

improves as e decreases, and exactly recovers the target posterior as e — 0, as then 
lim^^Q IT ^{y\x, 9) becomes a point mass at y = x (IReeves and Pettitt 20051) . 



Posterior simulation from TTLp{6\y) can then proceed via standard simulation al- 
gorithms, replacing pointwise evaluations of irLFi^^ly) with Monte Carlo estimates 
through the expectation (12.21) . based on draws x^, . . . ,x^ ~ 7r(a;|6') from the model 



(e.g. Marjoram et al. 2003 ). Alternatively, simulation from the joint posterior 7r^P'(6', x|?/) 
is available by contriving to cancel the intractable likelihoods 7r(x|6') in sample weights 
or acceptance probabilities. For example, importance sampling from the prior predic- 



tive distribution n^O, x) = 7r(x|6')7r(6') results in an importance weight of nLp{9, x\y)/Tx{6, x) oc 



'Ke{y\x^9), which is free of hkehhood terms. See Sisson et al. (2009) for a discussion 
of marginal and joint-space likelihood-free samplers. 

In general, the distribution of 'k{x\9) will be diffuse, unless x is discrete and dim(x) 
is small. Hence, generating x ~ 7r(-|^) with x ~ ?/ is improbable for realistic datasets 

and as a result the degree of computation required for a good likelihood-free ap- 
proximation TTLp{9\y) ^ 7r(^|y) (i.e. with small e) will be prohibitive. In practice, the 
function 7T^{y\x,9) is expressed through low dimensional vectors of summary statis- 
tics, >S'(-), such that 7r^{y\x,9) weights the intractable posterior through (12. ip with 
high values in regions where S{y) ~ S{x). 

If S{-) is sufficient for 9, then letting e — ^ recovers nLF{9\y) = TT{9\y) as before, 
but with more acceptable computational overheads, as dim(5'(x)) << dim(x). As suf- 
ficient summary statistics are generally unavailable, the use of non-sufficient statistics 
is commonplace. The effect of less efficient estimators of 9 in (12. 2p is a more diffuse 
approximation of 7r(6'|y). Hence the choice of summary statistics in any application 
is critical, with the ideal being low- dimensional, efficient and near-sufficient. 

In this article, we implement the likelihood-free sequential Monte Carlo sampler 



of Peters et al. (2008) , detailed in Appendix A. As the class of particle-based algo- 
rithms is the most efficient currently available in likelihood-free computation (e.g. 



McKinley et al. 2009 ), and within this class, the sampler of Peters et al. (2008) 



IS 



the only one to allow non-uniform functions 7r^{y\x, 9), this sampler provides the best 
combination of efficient simulation and flexible modelling. 
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3 Bayesian a-stable models 

We now develop univariate and multivariate Bayesian a-stable models. Unlike exist- 
ing methods, likelihood-free inference is independent of model parameterization. 



3.1 Univariate a-stable Models 

Denote the characteristic function of n i.i.d. univariate a-stable distributed random 
variables , . . . , X„ by $x (t) . A popular and convenient parameterization is 



$x it) 



exp - 7" [1 +2/?tan^sgn(t) (|7t|^"" - l)] ) if a ^ 1 
exp - 7 |t| [l + z/3fsgn(t)ln(7|t|)]) if a = 1, 



(3.1) 

where sgn (t) = ttt and = —1 (e.g. Samorodnitsky and Taqqu 1994 ). Many alter- 



native parameterizations are detailed in Nolan (2007) and Zolotarev (1986) Under 



(13. ip . the intractable stable density function is continuous and unimodal, taking sup- 
port on (— oo, 0) if a < 1, /3 = —1; (0, oo) if a < 1, = 1 and (— oo, oo) otherwise. 

Efficient simulation of auxiliary data, x ~ 7r{x\6), under the model is critical for 
the performance of likelihood-free methods (Section [2]). Here, it is straightforward 
to generate a-stable variates under the model defined by the characteristic function 
(13. ip (e.g. Devroye 1986: [Nolan 2007|) . This approach is provided in Appendix B. 



3.1.1 Summary statistics 

A key component of likelihood-free inference is the availability of low- dimensional, 
efficient and near-sufficient summary statistics. Since a-stable models can possess 
infinite variance [a > 1) and infinite mean (a < 1), this choice must be made with 
care. Here we present several candidate summary vectors, S1S5, previously utilized 
for parameter estimation in the univariate a-stable model. In Section H] we evaluate 
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the performance of these vectors, and provide informed recommendations for the 
choice of summary statistics under the hkehhood-free framework. 



Si McCuUoch's Quantiles 



McCuUoch (1986) and McCulloch (1998) estimate model parameters based on sample 
quantiles, while correcting for estimator skewness due to the evaluation of ^(x), the 

(0 



p*^ quantile of x, with a finite sample. Here, the data xu) are arranged in ascending 



order and matched with qs(i){x), where s (i) = Linear interpolation to p from 

the two adjacent s{i) values then establishes qp{x) as a consistent estimator of the 
true quantiles. Inversion of the functions 

^ _ 90.95 (•) - 90.05 (•) - _ 90.95 (•) +g0.05(-) - 2go.5(-) ^ _ go.75(-) - 90.25(-) 
qo.75{-) - qo.25{-) <?0.95(-) " go.OSl") 7 

then provides estimates of a, (3 and 7. Note that from a computational perspective, 
inversion of v^, vp or is not required under likelihood-free methods. Finally, we 
estimate 6 by the sample mean X (when a > 1). Hence Si{x) = {va,Vi3,v^,x). 

S2 Zolotarev's Transformation 



Based on a transformation of data from the a-stable family X — > Z, Zolotarev (1986) 
(p. 16) provides an alternative parameterization of the a-stable model (a,/3, 7,5) 
(z/, 7], t) with a characteristic function of the form 

log$z (i) = - exp |z/~^ log |t| + r - z^T^sgn (t) + C (^z/~^ - 1 j | , (3.2) 

where C is Euler's constant, and where i/ ^ i, \r]\ < min{l,2-yz/ — 1} and |r| < 
00. This parameterization has the advantage that logarithmic moments have simple 
expressions in terms of parameters to be estimated. For a fixed constant < < | 
(IZolotarev 19861 recommends ^ = 0.25) and for integer n/3, the transformation is 

Zj = Xsj-2 - e^3,-i - (1 - 0^3,, J = 1, 2, ... , n/3. 



Defining Vj = log \ Zj\ and Uj = sgn{Xj), estimates for i/, r/ and r are tlien given by 

i? = max{z>, (1 + |r/|)V4}, r/ = E (f/) , t = E{V), 

wfiere u = -^S^{V) — ^S^{U) + 1, using sample variances S'^{V) and S'^{U). As 
before, 6 is estimated by X (for a > 1), and so S2{x) = {d,rj,T,x). 

Ss Press's Method Of Moments 

For a ^ 1 and unique evaluation points ti,t2, ta, t^, the method of moments equations 
obtained from log$x(^) can be solved to obtain (IPress 1972( IWeron 2006P 

, log iog(- log mt2)\) - log iog(- log mt^)\) ^ logSS 



log \ti/t2 



U{ti) /ti-u{tz) /t^ - |t4r ^(tg) /ts - Ital" U{ti) /ti 

= 



^s,7 (|t,|-i-|t3|"-i)^3tan(f) 
where u(t) = tan~^[^"^]^ cos(tXj)/ ^"^j^ sin(txj)]. We adopt the evaluation points 



I, la— 1 I, iQ— 1 

r4| ~ Fsl 



ti = 0.2,^2 = 0.8,^3 = 0.1 and ^4 = 0.4 as recommended by Koutrouvelis (1980), and 
accordingly obtain 5'3(x) = (a, /3,7, 5). 

^4 Empirical Characteristic Function 

The empirical characteristic function, ^x{t) = ~ ^"=1 e**"^^ for t G (—00,00), can 
be used as the basis for summary statistics when standard statistics are not avail- 
able. E.g. this may occur through the non-existence of moment generating functions. 
Hence, we specify S4{x) = ($x(^i), • • • , *^'x(^2o)) where ti G {±0.5, ±1, ±1.5, . . . , ±5}. 

^5 Mean, Quantiles and Kolmogorov-Smirnov Statistic 

The Kolmogorov-Smirnov statistic is defined as KS{X) = sup^ \F^{z) — F^{z)\, the 
largest absolute deviation between the empirical cumulative distribution functions of 
auxiliary (X) and observed (Y) data, where F^{z) = ^ Y17=i hxi<z) and I{Xi<z) = 1 
ii Xi < z and otherwise. We specify S^i^x) = (x, {qp{x)} , K S (x)) , where the set of 
sample quantiles {qp{x)} is determined byp G {0.01, 0.05, 0.1, 0.15, ... , 0.9, 0.95, 0.99}. 
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Likelihood-free inference may be implemented under any parameterization which 
permits data generation under the model, and for which the summary statistics are 
well defined. From the above S1-S5 are jointly well defined for a > 1. Hence, to 
complete the specification of the univariate a-stable model we adopt the independent 
uniform priors a ~ U[l.l,2], p ~ U[-l,l], 7 ~ U[0,300] and 5 ~ U[-300,300] (e.g. 
IBuckle 19951) . Note that the prior for a has a restricted domain, reflecting the use of 
sample moments in and 5*5. For S4 we may adopt a ~ U(0, 2]. 



3.2 Multivariate a-stable Models 

Bayesian model specification and simulation in the multivariate a-stable setting is 
challenging (INolan et al. 200"T1 INolan 2008} Samorodnitsky and Taqqu 1994 ). Here 



we follow Nolan (2008) , who defines the multivariate model for the random vector 
X = {Xi, . . . , Xd) G M'^ through the functional equations 



-(t) = / |(t,s)rr(rfs); /5 (t) = a-" (t) / sgn((t,s))|(t,s)rr(c^s) 



(t,/x°) a^l 
(t,Ai°)-f 4(t,s)lnKt,s)|r(rfs) a = l 

where t = {ti, . . . ,td), s = (si, . . . , s^), {t,s) = J2i=i'^i^i^ ^'^ denotes the unit en- 
sphere, r (ds) denotes the unique spectral measure, and cr" (t) represents scale, P (t) 
skewness and /i (t) location (through the vector = {fi^, . . . , yuJJ)). Scaling properties 
of the functional equations, e.g. fi{rt) = r/i(t), mean that it is sufficient to consider 
them on the unit sphere (INolan 19971) . The corresponding characteristic function is 



$x (t) =Eexp (z (X,t)) =exp(-/x(t)+z(M°,t)) 
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with /x(t) = J^di^a ((t, s)) r((is), where the function is given by 
The spectral measure r(-) and location vector fj,^ uniquely characterize the mul- 



|n| 
Inl 



(l - isgn(u)tan (^)) a^l 
(l — z^sgn (m) In \u\) a = 1. 



tivariate distribution (Samorodnitsky and Taqqu 1994) and carry essential informa- 



tion relating to the dependence between the elements of X. The continuous spec- 
tral measure is typically well approximated by a discrete set of k Dirac masses 

k 



r (■) = ^WjSsji-) (e.g. [Byczkowski et al. 1993D where wj and 5sj(-) respectively 
denote the weight and Dirac mass of the j*'' spectral mass at location Sj G S"^. 
By simplifying the integral in /x(t), computation with the characteristic function 
= exp{— X]j=i Sj))} becomes tractable and data generation from the 

distribution defined by $x(^) efficient (Appendix B). As with the univariate case 
(13. ip . standard parameterizations of $x (t) ^i^^ be discontinuous at a = 1, result- 
ing in poor estimates of location and r(-). In the multivariate setting this is over- 
come by Zolotarev's M-parameterization (INolan et al. 200"T1 INolan 20080 . Although 
likelihood- free methods are parameterization independent, it is sensible to work with 
models with good likelihood properties. 

In a Bayesian framework we parameterize the model via the spectral mass, which 
involves estimation of the weights w = {wi, . . . ,Wk) and locations Si:fc G §'^^'^ of 
r(-). For k = 2 this corresponds to Sj = (cos</)j, sin^j) G S^. More generally, for 
Sj = (^s], . . . , sf) G B>'^, we use hyperspherical coordinates = {(f)\, . . . , where 

(pd-i = tan > • • • . 01 = tan M 

We define priors for the parameters (lu, 01:^, a), with </)i,k = {(f)i, . . . , 4>k), as 
w ~ Dirichlet(s, . . . , s), 0}, ~ U(0, 27r), /i° ~ N(^, a ~ U(0, 2) for z = 1, . . . , fc, 
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j = l,...,d, j' = l,...,d — l and impose the ordering constraint s}_^ < s] for all i on 
the first element of each vector Sj. Note that by treating the weights and locations 
of the spectral masses as unknown parameters, these may be identified with those 
regions of the spectral measure with significant posterior mass. This differs with the 



approach of Nolan (1997) where the spectral mass is evaluated at a large number 
of deterministic grid locations. In estimating r(-) | X a significant reduction in the 
number of required projection vectors is achieved (Sections 13.2.11 and 14. 2p . Further, 
the above prior specification does not penalize placement of spectral masses in close 
proximity. While this proved adequate for the presented analyses, alternative priors 



may usefully inhibit spectral masses at similar locations (IPievatolo and Green 1998]) . 



3.2.1 Summary statistics 

Sq Nolan, Panorska & McCulloch Projection Method 

For the d-vahaXe a-stable observations, Xj,2 = 1, . . . ,n, we take projections of Xj 
onto a unit hypersphere in the direction t G S'^. This produces a set of n univariate 
values, X* = {Xf, . . . ,X^), where = (Xj,t). The information in X* can then be 
summarized by any of the univariate summary statistics S'i(X*), . . . , 5*5 (X*). This 
process is repeated for multiple projections over ti, . . . , t,-. With the location param- 
eter estimated by x, for sufficient numbers of projection vectors, r, the summary 
statistics Sq(X) = (x, ^^(X*! ),..., ^^(X*-)) for some s e {1,2,3,4,5}, will capture 
much of the information contained in the multivariate data, if Sg is itself informative. 
The best choice of univariate summary vector Sg will be determined in Section 14.11 
We adopt a randomized approach to the selection of the projection vectors ti, . . . , t,-, 
avoiding curse of dimensionality issues as the dimension d increases f INolan 2008^ . 
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4 Evaluation of model and sampler performance 

We now analyze the performance of the Bayesian a-stable models and likelihood-free 
sampler in a sequence of simulation studies. For the univariate model, we evaluate 
the capability of the summary statistics Si, . . . , S5, and contrast the results with the 



samplers of Buckle (1995) and Lombardi (2007) The performance of the multivariate 



model under the statistics 5*6 is then considered for two and three dimensions. 

In the following, we implement the likelihood-free sequential Monte Carlo algo- 



rithm of Peters et al. (2008) (Appendix A) in order to simulate from the likelihood- 
free approximation to the true posterior 7rLp{6\y) ^ 7r(^|?/) given by (12. 2p . This 
algorithm samples directly from iTLFiGly) using density estimates based on P > 1 
Monte Carlo draws x^, . . . ~ 7r(x|6') from the model. We define 7r^{y\x,6) as a 
Gaussian kernel so that the summary statistics S{y) ~ N{S{x),e'^T,) for a suitably 
chosen S. All inferences are based on = 1000 particles drawn from 7rLF{6\y). Detail 
of algorithm implementation is removed to Appendix A for clarity of exposition. 



4.1 Univariate summary statistics and samplers 

We simulate n = 200 observations, y, from a univariate a-stable distribution with 
parameter values a = 1.7, (3 = 0.9, 7 = 10 and 6 = 10. We then implement the 
likelihood- free sampler targeting 7r^^(6'|i/) for each of the univariate summary statis- 
tics Si-S^ described in Section 13.1.11 with uniform priors for all parameters (Section 
13. ip . Alternative prior specifications were investigated (ILombardi 2007MNolan 199711 . 
with little impact on the results. 

Posterior minimum mean squared error (MMSE) estimates for each parameter, 
averaged over 10 sampler replicates are detailed in Table [H Monte Carlo standard 
errors are reported in parentheses. The results indicate that all summary vectors apart 



13 



from 5*5 estimate a and b parameters well, and for 7, 6*3 and 5*4 perform poorly. Only 
5*1 gives reasonable results for /9 and for all parameters jointly. Figure [H illustrates a 
progression of the MMSE estimates of each parameter using 6*1, from the likelihood- 
free SMC sampler output for each sampler replicate. As the sampler progresses, the 
scale parameter e decreases, and the MMSE estimates identify the true parameter 
values as the likelihood-free posterior approximation improves. 

The results in Table [1] are based on using P = 1 Monte Carlo draws from the 
model to estimate 'J^LF{Q\y) (cf. 12.21) within the likelihood-free sampler. Repeating 
the above study using P G {5, 10, 20} produced very similar results, and so we adopt 
P = 1 for the sequel as the most computationally efficient choice. 



For comparison, we also implement the auxiliary variable Gibbs sampler of Buckle (1995) 



and the MCMC inversion and series expansion sampler of Lombardi (2007) , based on 
chains of length 100,000 iterations (10,000 iterations burnin), and using their re- 
spective prior specifications. The Gibbs sampler performed poorly for most parame- 
ters. The MCMC method performed better, but has larger standard errors than the 
likelihood- free sampler using 5*1. 

The MCMC sampler ( ILombardi 20071) performs likelihood evaluations via inverse 
Fast Fourier transform (FFT) with approximate tail evaluation using Bergstrom ex- 
pansions. This approach is sensitive to a, which determines the threshold between 
the FFT and the series expansion. Further, as the tail becomes fatter, a finer spacing 
of FFT abscissae is required to control the bias introduced outside of the Bergstrom 
series expansion, significantly increasing computation. Overall, this sampler worked 
reasonably for a close to 2, though with deteriorating performance as a decreased. 
The Gibbs sampler (IBuckle 19951) performed extremely poorly for most settings and 
datasets, even when using their proposed change of variables transformations. As 
such, the results in Tabled] represent simulations under which both Gibbs and MCMC 



14 



samplers performed credibly, thereby typifying their best case scenario performance. 
4.2 Multivariate samplers 

We consider varying numbers of discrete spectral masses, k, in the approximation 

k 

to the spectral measure F (■) = Yl "^j^s {■)■ We assume that the number of spec- 

i=i 

tral masses is known a priori, and denote the d-variate a-stable distribution by 
Sa {d, k,w, <pi;k, n'^). Priors are specified in Section [321 Following the analysis of 
Section l4m we incorporate 5*1 within the summary vector Sq. 

For datasets of size n = 200, we initially consider the performance of the bivariate 
a-stable model, Sa (2, k, w, (pi-.k, 0), for /c = 2 and 3 spectral masses, with respect to 
parameter estimation and the impact of the number of projection vectors ti, . . . , t,-. 

The true and mean MMSE estimates of each parameter, placing projection vectors 
at the true locations of the spectral masses, are presented in Table O In addition, 
results are detailed using t = 2,5, 10 and 20 randomly (uniformly) placed projection 
vectors, in order to evaluate the impact of spectral mass location uncertainty. The 
likelihood-free sampler output results in good MMSE parameter estimates, even for 
2 randomly placed projection vectors. The parameter least accurately estimated is 
the location vector, Directly summarized by a sample mean in Sq{-), estimation 
of location requires a large number of observations when the data have heavy tails. 

Figure [2] illustrates progressive sampler performance for the Sa (2, 2, w, (pi:2, 0) 
model, with a = 1.7, w = (7r/4,7r) and (pi-^ = (7r/4,7r). Each circular scatter plot 
presents MMSE estimates of weight (radius) and angles (angle) of the two spectral 
masses, based on 10 sampler replicates. The sequence of plots (a)-(d) illustrates 
the progression of the parameter estimates as the scale parameter e (of 'K^{y\x,6)) 
decreases (and hence the accuracy of the likelihood- free approximation 7rLF{0\y) ~ 
Tc{6\y) improves). As e decreases, there is a clear movement of the MMSE estimates 
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towards the true angles and weights, indicating appropriate sampler performance. 

With simulated datasets of size n = 400, we extend the previous bivariate study 
to 3 dimensions, with k = 2 discrete spectral masses. The true parameter values, 
and posterior mean MMSE estimates and associated standard errors, based on 10 
sampler replicates, are presented in Table [3l Again, reasonable parameter estimates 
are obtained (given finite data), with location again the most difficult to estimate. 

In analogy with Figure [21 progressive sample performance for the first spectral 
mass (with wi = 0.7 and (pi = (7r/4, tt)) for decreasing scale parameter e is shown 
in Figure O Based on 200 replicate MMSE estimates (for visualization purposes), 
the shading of each point indicates the value of Wi as a percentage (black=0%, 
white=100%), and the location on the sphere represents the angles (pi. For large 
e, the MMSE estimates for location are uniformly distributed over the sphere, and 
the associated weight takes the full range of possible values, — 100%. As e decreases, 
the estimates of spectral mass location and weight become strongly localized and cen- 
tered on the true parameter values, again indicating appropriate sampler performance. 
Similar images are produced for the second discrete spectral mass. 

5 Analysis of exchange rate daily returns 

Our data consist of daily exchange rates for 5 different currencies recorded in GBP 
between 1 January 2005 and 1 December 2007. The data involve 1065 daily-averaged 
LIBOR (London interbank offered rate) observations y[, . . . , ^loes- '^^^ standard log- 
transform generates a log returns series yt = In {y[j^i/y'f) . Cursory examination of each 
returns series reveals clear non-Gaussian tails and/or skewness (Table HI bottom). 

We initially model each currency series as independent draws from a univariate a- 
stable distribution. Posterior MMSE parameter estimates for each currency are given 
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in Table m based on 10 replicate likelihood-free samplers using the 5*1 summary vec- 
tor. For comparison, we also compute McCuUoch's sample quantile based estimates 
(derived from Si, c.f. Section 13.1. ip . and maximum likelihood estimates using J. 
P. Nolan's STABLE pTogram (available online), using the direct search SPDF option 
with search domains given by a G (0.4, 2], (3 e [-1, 1], 7 G [0.00001, 1] and 6 G [-1, 1]. 
Overall, there is good agreement between Bayesian, likelihood- and sample-based es- 
timators. All currency returns distributions are significantly different from Gaussian 
(a = 2,/3 = 0), and exhibit similar family parameter (a) estimates over this time pe- 
riod. However, the GBP to YEN conversion demonstrates a significantly asymmetry 
{(3) compared to the other currencies. 

An interesting difference between the methods of estimation, is that McCuUoch's 
estimates of a differ considerably from the posterior MMSE estimates, even though 
the latter are constructed using McCuUoch's estimates directly as summary statistics. 
Si. One reason that the Bayesian estimates are more in line with the MLE's, is that 
likelihood-free methods largely ignore bias in estimators used as summary statistics 
(comparing the closeness between biased or unbiased estimators will produce similar 
results - consider comparing sample and maximum likelihood estimators of variance). 

The multivariate a-stable distribution assumes that its marginal distributions, 
which are also a-stable, possess identical shape parameters. This property implies 
important practical limitations, one of which is that it is only sensible to jointly model 
data with similar marginal shape parameters. Accordingly, based on Table HI we now 
consider a bivariate analysis of AUD and EURO currencies. Restricting the analysis to 
the bivariate setting also permits comparison with the bivariate frequentist approach 



described in Nolan (1997) using the MVSTABLE software (available online). 

A summary illustration of the discrete approximations to the underlying continu- 
ous spectral mass is shown in Figure HI Assuming /c = 3 discrete spectral masses and 
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based on 10 likelihood-free sampler replicates, the mean MMSE posterior estimates 
(solid black line) with mean So" posterior credibility intervals (dotted line), identify re- 
gions of high spectral mass located at 2.7, 3.9 and 5.6, with respective weights 0.45, 0.2 



and 0.35. Broken lines in Figure H] denote the frequentist estimates of Nolan (1997) 



based on the identification of mass over an exhaustive mesh grid using 40 (dashed 
line) and 80 (dash-dot line) prespecified grid locations (projections). 

Overall, both approaches produce comparable summary estimates of the spectral 
mass approximation, although the likelihood-free models generate full posterior distri- 
butions, compared to Nolan's frequentist estimates. The assumption of A; = 3 discrete 
spectral masses provides a parsimonious representation of the actual spectral mass. 
For example, the spectral mass located at 2.7 accounts for the first two/three masses 
based on Nolan's estimates (80/40 projections). While the frequentist approach is 
computationally restricted to bivariate inference, the likelihood-free approach may 
naturally be applied in much higher dimensions. 



6 Discussion 

Statistical inference for a-stable models is challenging due to the computational in- 
tractability of the density function. In practice this limits the range of models fitted, 
to univariate and bivariate cases. By adopting likelihood-free Bayesian methods we 
are able to circumvent this difficulty, and provide approximate, but credible posterior 
inference in the general multivariate case, at a moderate computational cost. Critical 
to this approach is the availability of informative summary statistics for the parame- 
ters. We have shown that multivariate projections of data onto the unit hypersphere, 
in combination with sample quantile estimators, are adequate for this task. 

Overall, our approach has a number of advantages over existing methods. There is 
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far greater sampler consistency than alternative samplers, such as the auxiliary Gibbs 
or MCMC inversion plus series expansion samplers (IBuckle 1995| ILombardi 20071) . It 



is largely independent of the complexities of the various parameterizations of the a- 
stable characteristic function. The likelihood-free approach is conceptually straight- 
forward, and scales simply and is easily implemented in higher dimensions (at a higher 
computational cost). Lastly, by permitting a full Bayesian multivariate analysis, the 
component locations and weights of a discrete approximation to the underlying con- 
tinuous spectral density are allowed to identify those regions with highest posterior 
density in a parsimonious manner. This is a considerable advantage over highly com- 
putational frequentist approaches, which require explicit calculation of the spectral 
mass over a deterministic and exhaustive grid (e.g. INolan 19971) . 



Each analysis in this article used many millions of data-generations from the 
model. While computation for likelihood-free methods increases with model dimen- 
sion and desired accuracy of the model approximation (through e), much of this can 
be offset through parallelization of the likelihood- free sampler (IPeters et al. 20081) . 

Finally, while we have largely focused on fitting a-stable models in the likelihood- 
free framework, extensions to model selection through Bayes factors or model averag- 
ing are immediate. One obvious candidate in this setting is the unknown number of 
discrete spectral masses, /c, in the approximation to the continuous spectral density. 
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Appendix A 

SMC sampler PRC-ABC algorithm (IPeters et al. 2008|) 



Initialization: Set t = 1 and specify tolerance schedule ei, . . . , er- 

For z = 1, . . . , iV, sample ~ 7r(^), and set weights W^{ef) = 7rLF,i{d?\y)/Tr{e?) 

Resample: Resample TV particles with respect to WtiOl^'') and set VF'f(^f^) = j^, 

2 = 1,. ..,iV. 

Mutation and correction: Set t = t + 1 and i = 1: 

(a) Sample ~ Mt{9t) and set weight for of^ to 

Wt{ef^) = nLF.M\y)/Mt{e?). 

(b) With probability 1 — p^*^ = 1 — min{l, Wt{6f^)/ct}, reject 9^^ and go to (a) 

(c) Otherwise, accept Of^ and set Wt{ef^) = Wt{eP)/p^'l 

(d) Increment i = i + 1. If z < A^, go to (a). 

(e) If t < T then go to Resample. 

This algorithm samples weighted particles from a sequence of distributions nLF,t{^\y) 
given by (12.21) . where t indexes a sequence of scale parameters ei > ... > e^. 
The final particles {{Wt{o!^^),oP) : i = 1, ...,A^}, form a weighted sample from 
the target vrLF,T(^|y) (e.g. IPeters et al. 20081) . The densities hlfA^Iv) are esti- 
mated through the Monte Carlo estimate of the expectation (12. 2p based on P draws 
x^, . . . , ~ 7r(x|6'). 

For the simulations presented we implement the following specifications: for uni- 
variate a-stable models 6 = [a, /5, 7, 6) and for multivariate models 6 = [w, cf>i:k, M*^, a); 
we use = 1000 particles, initialized with samples from the prior; the function 
7r^(y\x,9) is defined by S{y) ~ N(S{x), e'^T.) where E is an estimate of Cov(S'(x)|^) 
based on 1000 draws x^, . . . , ~ 7r(x|6') given an approximate maximum like- 



lihood estimate 6 of 6 ( Jiang and TurnbuU 2004 ); the mutation kernel Mt{6t) 
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wl;l!i{6j:^}i)(f){6t] Ol^li, A) is a density estimate of the previous particle population 
{(W^t-i (^1^1)7 ^l-i) : « = 1, • • • , N}, with a Gaussian kernel density with covariance 
A; for univariate a-stable models A = diag(0.25, 0.25, 1,1), and for multivariate mod- 
els A = diag(l, . . . , 1, 0.25) (with Dirichlet proposals and kernel density substituted 
for w); the sampler particle rejection threshold is adaptively determined as the 90*^* 
quantile of the weights Cf = qo.9{{wl;^\6l^^)}) , where {wl;^\6l^^)} are the particle 
weights prior to particle rejection (steps (b) and (c)) at each sampler stage t (see 
IPeters et al. 20081) . 

For each analysis we implement 10 independent samplers (in order to monitor algo- 
rithm performance and Monte Carlo variability), each with the deterministic scale pa- 
rameter sequence: G {1000, 900, . . . , 200, 100, 99, . . . , 11, 10, 9.5, 9, . . . , 5.5, 5, 4.95, 
. . . , 3.05, 3, 2.99, 2.98, . . . , 0.01, 0}. However, we adaptively terminate all samplers at 
the largest e value such that the effective sample size (estimated by [Xlili [^t^ (^f ^)]^] ""^y 
consistently drops below 0.2A^ over all replicate sampler implementations. 

Appendix B: Data generation 



Simulation of univariate a-stable data {DuMouchel 1975, Chambers et al. 1976 ) 

1. Sample W ~ Exp(l) to obtain w 

2. Sample U ~ Uniform[— 7r/2, 7r/2] to obtain u 

3. Apply transformation to obtain sample y 



^a,f3 (cos w)"/2 



cos(«— a('u+ijQ,/3)) 



(1 + Pu) tanu - /51n ^T^"^^" 



if a ^ 1 
if a = 1 



with S^^f^ = (1 + /^^tan^ (^)) ^^^^ and B^^p = ^ arctan {(3tan (^)) . In this 
case y will have distribution defined by $x (t) with parameters {a, (3,1,0). 
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4. Apply transformation to obtain sample y = + S with parameters (a, (3, 7, S). 



Simulation of d- dimensional, multivariate a-stable data (Nolan 2001) 



1. Generate Zi, Zk i.i.d. random variables from the univariate a-stable distri- 
bution with parameters (a, /3, 7, 5) = (a, 1, 1, 0). 

2. Apply the transformation 

E ^/"Z,s, + /xO a ^ 1 

i=i 

E^.(^. + ^lnK))s,+/x° a = l 



with Si, . . . , Sfc, fi^ G S*^. Note that while the complexity for generating realizations 
from a multivariate a-stable distribution is linear in the number of point masses (k) 
in the spectral representation per realization, this method is strictly only exact for 
discrete spectral measures. 
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Lombardi 
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55 


a (1.7) 


1.77 (0.18) 


1.62 (0.10) 


1.69 (0.06) 


1.65 (0.07) 


1.70 (0.06) 


1.71 (0.04) 


1.56 (0.05) 


/3 (0.9) 


0.54 (0.21) 


0.86 (0.18) 


0.86 (0.10) 


0.65 (0.13) 


0.31 (0.09) 


0.38 (0.12) 


0.49 (0.11) 


7 (10.0) 


18.17 (6.19) 


9.59 (2.16) 


9.79 (0.21) 


10.44(0.56) 


38.89 (6.34) 


39.12 (5.92) 


9.34 (0.14) 


5 (10.0) 


12.30 (4.12) 


9.70 (2.19) 


10.64 (0.83) 


9.31 (0.86) 


10.25 (0.98) 


10.83 (1.34) 


11.18 (1.05) 



Table 1: Means and standard errors (in parentheses) of posterior MMSE estimates of a, 
/3, 7 and 5 under the univariate a-stable model, based on 10 sampler replicates. Parameter 
values used for data simulation are given in the left column. Comparisons are between the 
auxiliary variable Gibbs sampler method of Buckle (1995), the inversion MCMC method of 
Lombardi (2007), and the likelihood-free method, using summary statistics SiS^. 
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Projection vectors at locations of true spectral masses. 
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1.66 (0.04) 


0.16 (0.19) 


0.81 (0.65) 3.19 (0.46) 


0.55 (0.06) 


0.45 (0.05) 
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1.79 (0.02) 


0.36 (0.18) 


0.84 (0.27) 3.18 (0.29) 4.91 (0.24) 


0.35 (0.05) 


0.25 (0.04) 


0.40 (0.05) 
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1.67 (0.06) 


-0.13 (0.16) 


0.73 (0.55) 3.58 (0.57) 


0.58 (0.09) 


0.42 (0.10) 
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1.76 (0.05) 


-0.16 (0.26) 


0.91 (0.66) 3.65 (0.62) 4.85 (0.55) 


0.36 (0.10) 


0.24 (0.09) 


0.40 (0.08) 








5 projection vectors 
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1.71 (0.05) 


0.08 (0.14) 


0.71 (0.60) 3.80 (0.67) 


0.60 (0.07) 


0.40 (0.09) 
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1.75 (0.04) 


0.29 (0.17) 


0.86 (0.62) 3.68 (0.52) 4.82 (0.41) 


0.35 (0.09) 


0.20 (0.07) 


0.42 (0.09) 








10 projection vectors 
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1.72 (0.02) 


0.21 (0.21) 


0.75 (0.32) 3.31 (0.32) 


0.59 (0.05) 


0.41 (0.07) 
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1.73 (0.03) 


0.25 (0.16) 


0.76 (0.44) 3.31 (0.48) 4.78 (0.19) 


0.34 (0.07) 


0.24 (0.04) 


0.42 (0.05) 








20 projection vectors 
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1.71 (0.03) 


-0.14 (0.14) 


0.76 (0.36) 3.21 (0.23) 


0.63 (0.03) 


0.37 (0.03) 
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1.72 (0.02) 


0.18 (0.23) 


0.77 (0.32) 3.25 (0.31) 4.75 (0.15) 


0.34 (0.04) 


0.23 (0.03) 


0.43 (0.03) 



Table 2: Mean MMSE parameter estimates (and standard errors) for the bivariate 
a-stable Sa (2, k, w, <pi;k, fi^) model, for k = 2,3 discrete spectral masses, calculated 
over 10 replicate samplers. Projections vectors are placed at the true, and 2, 5, 10 
and 20 randomly selected spectral mass locations. The true value of is the origin. 
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2 1.71 (0.02) 


0.53 (0.89) 


1.12 (0.34) 


3.81 (0.45) 1. 


34 (0.54) 


4.24 (0.69) 


0.28 (0.06) 


0.72 (0.05) 



Table 3: Mean MMSE parameter estimates (and standard errors) for the trivariate 
a-stable Sa (3, 2, w, (pi,2, /x°) model, with k = 2 discrete spectral masses, calculated 
over 10 replicate samplers. The true value of fi^ is the origin. 
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1.39 


1.38 
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quantilc 
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Skcwness 
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-0.03 
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Std. dev. 




0.004 
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0.004 


0.003 


Mean 




-4e-5 


-4e-5 
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Table 4: Posterior MMSE estimates (Monte Carlo errors) from the likelihood-free 
model, and maximum likelihood estimates (standard deviation). MLE's, parameter 
estimates using McCuUoch's quantile (McCuUoch, 1998), and sample statistics (mean, 
standard deviation, skewness and kurtosis) obtained from J. P. Nolan's STABLE 
software, available at academic2.american.edu/~jpnolan. 
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Figure 1: Traces of posterior MMSE estimates of a, /J, 7 and 6 under the univariate 
a-stable model and likelihood- free sampler (summary statistics based on 10 sampler 
replicates. Traces are shown as a function (x-axis) of sampler progression (t) and scale 
parameter reduction et < et-i- Parameter values used for data generation are a = 1.7, 
P = 0.9, 7 = 10 and 6 = 10. 
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Figure 2: Circular scatter plot of MMSE estimates for k = 2 spectral mass angles 
(angle) and weights (radius) for bivariate ct-stable Sa (2, 2, w, 4>i-2i 0) model, with 
a = 1.7, w = (0.4,0.6) and (pi-2 = (7r/4,7r). Plots (a)-(d) demonstrate evolution of 
the estimates for decreasing scale parameter values e, based on 10 sampler replicates. 
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Figure 3: Spherical heat map of MMSE estimates for the first of /c = 2 discrete 
spectral masses, 0i, for the trivariate a-stable Sa (3, 2, w, (f>i;2, 0) model. True values 
of the first spectral mass are Wi = 0.7 (70%) and (pi = (7r/4, tt). Point shading 
indicates MMSE value of Wi as a percentage. The plots demonstrate the evolution of 
the estimates for decreasing scale parameter values e, based on 200 sampler replicates. 



32 



■ ■ ■ Nolan (40 projections] 



■ I ■ I Nolan (80 projections) 




- 



J \ \ 

1 3 3 

Location of spectral mass 



Figure 4: Estimates of spectral mass location (x-axis) and cumulative weight (y- 
axis) for AUD and EURO currencies data. Solid line denotes mean posterior MMSE 
estimates of likelihood- free SMC sampler output, and dotted line illustrates mean 
3cr posterior credibility intervals, based on A; = 3 discrete spectral masses, 20 ran- 
domly placed projection vectors and 10 replicate samplers. Broken lines denote 
estimate of spectral mass using J. P. Nolan's MVSTABLE software, available at 
academic2 . american . edu/~ jpnolan, with (dashed line) 40 deterministic pro- 
jection locations and (dash-dot line) 80 projection locations. 
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